On the Eigenvalue Density of Real and Complex Wishart Correlation Matrices 
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Wishart correlation matrices are the standard model for the statistical analysis of time series. 
The ensemble averaged eigenvalue density is of considerable practical and theoretical interest. For 
complex time series and correlation matrices, the eigenvalue density is known exactly. In the real 
case, however, a fundamental mathematical obstacle made it forbidingly complicated to obtain exact 
results. We use the supersymmetry method to fully circumvent this problem. We present an exact 
formula for the eigenvalue density in the real case in terms of twofold integrals and finite sums. 

PACS numbers: 05.45.Tp,02.50.-r,02.20.-a 

Keywords: Wishart correlation matrices, eigenvalue density, supersymmetry 



Time series analysis is an indispensable tool in the 
study of complex systems with numerous applications in 
physics, climate research, medicine, signal transmission, 
finance and many other fields [H-Q. A time series is an 
observable such as the water level of a river, the tem- 
perature, the intensity of transmitted radiation, the neu- 
ron activity in electroencephalography (EEG), the price 
of a stock, etc., measured at usually equidistant times 
k = 1,. ..,n. Suppose we measure p such time series 
Mj, j = 1, . . . ,p, for example, in the case of EEG, at p 
electrodes placed on the scalp or, in the case of tempera- 
tures, at p different locations. Our data set then consists 
of the pxn rectangular matrix M with entries M.-fc. The 
time series Mj are usually real (labeled (3 = 1), but in 
some applications they can be complex (f3 — 2). Often, 
one is interested in the correlations between the time se- 
ries. To estimate them, one normalizes the time series 
Mj to zero mean and unit variance. The correlation co- 
efficient between the time series Mj and Mi is then given 
as the sample average 
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= - J2 M ikMf k and 



k=l 



C = -MM^ (1) 
n 



is the correlation matrix. For real time series (j3 = 1), 
the complex conjugation is not needed and the adjoint 
is simply the transpose. We notice that C is a p x p 
real-symmetric (j3 = 1) or Hcrmitean (Jj — 2) matrix. 

The eigenvalues of C provide important information, 
see recent examples in Refs. [4, 5]. As the empirical infor- 
mation is limited, it is desirable to compare the measured 
eigenvalue density with a "null hypothesis" that results 
from a statistical ensemble. The ensemble is defined |6j 
by synthetic real or complex time series Wj, j — 1, . . . ,p 
which yield the empirical correlation matrix C upon av- 
eraging over the probability density function 



P^(W,C)~exp 



w 



that is, we have by construction 

I d\W]Pa{W,C)-WW"' = C 
J n 
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where the measure d[W] is the product of the differentials 
of all independent elements in W. To ensure that C is in- 
vertible, we always assume n > p. When going to higher 
order statistics, the Gaussian assumption ([2]) is not nec- 
essarily justified, but it often is a good approximation. 
This multivariate statistical approach is closely related 
to Random Matrix Theory and the matrices WW^ 
are referred to as Wishart correlation matrices. One is 
interested in the ensemble averaged eigenvalue density of 
these matrices. In terms of the resolvent, it reads 

S p {x) = -±*nf d[Wn(W,C)tr x+ip \ wV (4) 

where t p is the pxp unit matrix. The argument x carries 
a small positive imaginary increment e > 0, indicated by 
the notation x + = x + ie. The limit s — > is implicit 
in our notation. Due to the invariance of the trace and 
the measure, the ensemble averaged eigenvalue density 
Sp(x) only depends on the eigenvalues Aj, j — 1, . . . ,p 
of C. Hence we may replace C in Eq. (j4|) by the pxp 
diagonal matrix A = diag(Ai, . . . , A p ). We notice that 
the eigenvalues are positive definite, Aj > 0. 

A large body of literature is devoted to the eigenvalue 
density (j4|). Its asymptotic form for large n and p has 
been studied in great detail, see Refs. [8, 9]. However, an 
exact closed-form result for finite n and p is only avail- 
able in the complex case 2, 11 1. Unfortunately, a deep, 
structural mathematical reason made it up to now impos- 
sible to derive such a closed-form result in the real case 
which is the more relevant one for applications. We have 
three goals: We, first, introduce the powerful supersym- 
metry method 12|-|l4l] to Wishart correlation matrices for 
arbitrary C. This has, to the best of our knowledge, not 
been done before. We, second, use the thereby achieved 
unique structural clearness to derive a new and exact 
closed-form result for the eigenvalue density in the real 
case for finite n and p. We, third, show that our results 
are easily numerically tractable and compare them with 
Monte Carlo simulations. 

Why does the real case pose such a substantial prob- 
lem? — This is best seen by going to the polar decom- 
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position W = UwV, where U G 0(p), V G 0(n) for 
/3 = 1 and U G U(p), G U(n) for ,8 = 2 and where 
it) is the p x n matrix containing the singular values or 
radial coordinates Wj, j = I,..., p. In particular, one 



has WW^ = Uw 2 W, with w 2 = 



When inserting 



into Q, one sees that the non-trival group integral 
$^(A,w; 2 ) = j exp (-^trU^A^UwA d/j,(U) 



(5) 



has to be done to obtain the joint probability density 
function of the radial coordinates Wj. Here, dp,(U) is 
the invariant Haar measure. For ft = 2, this integral 
is the celebrated Harish-Chandra-Itzykson-Zuber inte- 
gral and known explicitly 7jl 16 1- For ft = 1, however, 
4> 1 (A,u' 2 ), is not a Harish-Chandra spherical function, it 
rather belongs to the Gelfand class 17[ and a closed-form 
expression is lacking. The only explicit form known is a 
cumbersome, multiple infinite series expansion in terms 
of zonal or Jack polynomials @, [l8| . This inconvenient 
feature then carries over to the eigenvalue density (j4j, 
but we will arrive at a finite series over twofold integrals. 
The supersymmetry method is based on writing 



Sp{x) 



.j_ Im 9MJ) 

up oJ 



(6) 



j=o 



as the derivative of the generating function 



Z P {J) = / d[W\P p (W,K) 



det(x+t p + Jt p -WW^) 
det(x+l p - WW^) 



(7) 



with respect to the source variable J at J = 0. One has 
the normaliziation Zp(Q) = 1 at J = 0. We consider 
the real and the complex case and use the latter as test. 
We map Zp{J) onto superspace using steps which are by 
now standard, see Refs. TJl 14 1. A particularly handy 



approach for applications such as the present one is given 
in Ref. we use the same conventions and find 

Zp(J) = [ d[p]Ip{p) 



: ]~J sdct 



-/3/2 



"I4//3 - Jl ~ 2 Aj75 



(8) 



The merit of this transformation is the drastic reduc- 
tion in the number of degrees of freedom, because the 
variables to be intergrated over form the 4//3 x 4//3 Her- 
mitean supermatrix 



P11 
P12 



Pl2 
IP22 



(9) 



For /3 = 2, pn and p 22 are scalar, real commuting vari- 
ables and P12 is a complex anticommuting scalar variable. 
For /3 = 1, pn is a 2 x 2 real symmetric matrix, p 22 has 
to be multiplied with 1 2 and we have 



P12 



X X 



(10) 



where x,£ and x*i£* denote anticommuting variables 
and their complex conjugates, respectively. We also in- 
troduced the matrix 7 = diag(0 2 //3, — 12//3) and the su- 
persymmetric Ingham-Siegel integral 

Ip{p) = J d[o-]sdet~ n,3/2 (l 4/( 3 + ict) cxp(istrcrp) , (11) 

where a has the same form as p. The supertrace and 
superdeterminant [20I ] are denoted by str and sdet. 

Starting from the generating function (J8J) we first con- 
sider the complex case /3 = 2. By introducing eigenvalue- 
angle coordinates for the supermatrix p, we rederive in a 
straightforward calculation the eigenvalue density £2(2;) 
as found in Ref. [1(3]. In the real case j3 = 1, the 
analogous approach leads to inconvenient Efetov-Wegner 
terms 14j , and we thus proceed differently. Since the gen- 
erating function remains invariant under rotations of the 
matrix pn, we introduce its eigenvalues R\ — diag(ri, r 2 ) 
and the diagonalizing angle as new coordiantes. This 
yields the Jacobian |A 2 (i?i)| = \r\ — r 2 |. The next step 
is to evaluate the Ingham-Siegel integral Ii(p)- The su- 
permatrix a in Eq. (jlip has the same form as p in Eq. (0) . 
Doing the integral over a followed by an expansion in the 
anticommuting variables of p according to Eq. (|10j) gives 



h (p) ~ det ( "- 1)/2 i?! exp (-strp) 6(i?i) 

n— 2 / * \ / o \ n — 1 




xT + ZC 

n r 2 



i<9p : 
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idp : 



22 



U(P22) • 



(12) 



In a simple, direct calculation, we also expand the prod- 
uct of the superdeterminats in Eq. (jHJ in the anticom- 
muting variables of p. We collect everything and do the 
integration over the anticommuting variables. With the 
notation Qj = x + — 2Ajip 2 2, we obtain 

z x {j) ~ y y d P22 |A 2 (i? 1 )idet("- i )/ 2 i? 1 e(i? 1 ) 

p ( J + I 

exp(-(n + r 2 -2^ 2 ))n detl/2(a;+ ^ 

f Aet~ x Rx 



(J + Qj) 

, , „. , i (x+l 2 -2A j r 1 ) 



idp 



<) + £ (2A,) 2 / I 



t (J + Qj) V> + -2AjTi)r 2 



(2;+ — 2Ajr 2 )ri / \idp 
p 

E 



j> '2.)' 1 / \ "^H22 
t0 A j 



(2A^ 2 



(J + Q,)(.t+ -2A i n) 



(2A fc ) 2 



U(p2 2 ) -(13) 



(J + Q fc )(x+ -2A fe r 2 ) V^p 22/ 
According to Eq. ^ we have to take the derivative with 
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FIG. 1: Eigenvalue density for p = 5 and n = 200: analytical formula (solid lines) and Monte Carlo simulations (histogram 
with bin width 3). 



respect to J. This leads to the three relations 







v 



dJ u-l 

d 1 



= Ep-i{Q) 



(14) 



t[(J + Qi) 



J=0 



<9 1 1 ^ 



- -Sp-2y (0) 



Z=l 



,7=0 



where -B m -i,fc(Q) denotes the elementary symmetric poly- 
nomial of order m in the variables Qj, j = 1, . . . ,p ^ i,k 
with and Qfc omitted, 



E m ;i,k(Q) = 

We finally arrive at 



(15) 



l<-ii <-<i m <p 



S 1 ( X )=clraJd[R l] Jdp 22 uUdetl/2{x+l 



6(i?i)e- (ri+r2 ~ 2lp22) det -1 .Ri£p_i(Q) 



2 - 
9 



-Ew 2 



i 



(a; + - 2A 3 T2)ri 
1 

(x+ - 2A i n)r 2 



(16) 



E 



d 



(2A J ) 2 (2A fc ) 2 £; p _3, J , fc (Q) 
(a:+ - 2AjTi)(aJ+ - 2A fe r 2 ) Vidp 



Z<9/9 2 
n— 2\ 



where the constant reads c = (— l) n+1 / (Anp(n — 2)!). 
Due to the ^-distribution, the integral over p 22 are el- 
ementary. Hence we end up with an expression for the 
eigenvalue density S\{x) which essentially is a twofold 
integral. 



The integrals in Eq. (|16p can be numerically evaluated 
by using a regularisation technique of the type 



oo oo 



Im / dr\ I dr 2 



f(n,r 2 ) 



i=i 



!=1 v 
Ar 1 A" 1 



E 



dri / dr 2 f(n,r 2 ) 



0<i,j<P + .-i 
(i+j)G2N+l A i+l A j + l 



1 



n v^-ArV^-V 



(17) 



Here we assume an ordering of the eigenvalues such that 
A > Ai > . . . > A p > A p+ i with A 1 = and 
Al^j = oo. The real function /(ri,r 2 ) is independent 
of e and has no singularities. The singularties at the 
boundaries of the domain are integrable. Using the com- 
mercial software Mathematica® [21|, we evaluate our 
formula (1171) numerically. For independent comparison, 
we also carry out Monte Carlo simulations with ensem- 
bles of 10 5 matrices. In Figs. [1] and [2j we show the 
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results for p = 5 and p = 10 and n = 200 with 
the chosen empirical eigenvalues Aj, j = 1,...,5 
of {1.44,0.64,0.49,0.25,0.16} and A,-, j = 1, . . . , 10 of 
{1, 0.81, 0.7225, 0.64, 0.45, 0.36, 0.25, 0.2025, 0.1225, 0.03}, 
respectively. The agreement is perfect. 

In conclusion, we introduced the supersymmetry 
method for the first time to Wishart correlation matrices. 
We thereby derived exact expressions for the eigenvalue 
density in terms of low-dimensional integrals. This is a 
drastic reduction, as the original order of integrals is np. 
Our approach solves a serious mathematical obstacle in 
the real case. A presentation for a mathematics audience 



will be given elsewhere [22|. Here, we derived and dis- 
cussed the formulae needed for applications. In the real 
case (/3 = 1), we obtained the previously unknown exact 
solution in terms of a finite sum of twofold integrals. We 
evaluated our formula numerically and confirmed it by 
comparing to Monte Carlo simulations. 
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